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1.  Introduction 


We  are  concerned  in  this  report  with  tracking  moving  sources  using  a  net¬ 
work  of  aeroacoustic  sensors.  We  assume  that  the  sensors  are  placed  in  an 
"array  of  arrays"  configuration  containing  several  small-aperture  arrays 
distributed  over  a  wide  area.  Each  array  contains  local  processing  capa¬ 
bility  and  a  communication  link  with  a  fusion  center.  A  standard  approach 
for  estimating  the  source  locations  involves  bearing  estimation  at  the  in¬ 
dividual  arrays,  communication  of  the  bearings  to  the  fusion  center,  and 
processing  of  the  bearing  estimates  at  the  fusion  center  with  a  tracking  al¬ 
gorithm  [1-5].  This  approach  is  characterized  by  low  communication  band¬ 
width  and  low  complexity,  but  the  localization  accuracy  may  be  inferior  to 
the  optimal  solution  in  which  the  fusion  center  jointly  processes  all  of  the 
sensor  data.  The  optimal  solution  requires  high  communication  bandwidth 
and  high  processing  complexity.  The  amount  of  improvement  in  localiza¬ 
tion  accuracy  that  is  enabled  by  greater  communication  bandwidth  and 
processing  complexity  is  dependent  on  the  scenario,  which  we  character¬ 
ize  in  terms  of  the  source  motion  parameters,  the  power  spectra  (and  band¬ 
width)  of  the  signals  and  noise  in  the  sensor  data,  the  coherence  between 
the  source  signals  received  at  widely  separated  sensors,  and  the  observa¬ 
tion  time  (amount  of  data).  We  present  a  framework  to  identify  scenarios 
that  have  the  potential  for  improved  localization  accuracy  relative  to  the 
standard  bearings-only  tracking  method.  We  propose  an  algorithm  that  is 
bandwidth-efficient  and  nearly  optimal  that  uses  beamforming  at  small- 
aperture  sensor  arrays  and  time-delay  estimation  (TDE)  between  widely 
separated  sensors. 

The  sensor  signals  are  modeled  as  Gaussian  random  processes,  allowing 
the  inclusion  of  deterministic  as  well  as  random  propagation  effects.  We 
have  previously  studied  this  model  [6]  for  the  case  of  a  single  source  with 
fixed  position  (no  motion).  We  extend  the  analysis  in  this  report  to  moving 
sources  that  follow  a  parametric  motion  model. 

Accurate  time-delay  (TD)  estimates  are  required  in  order  to  achieve  im¬ 
proved  localization  accuracy  relative  to  bearings-only  triangulation.  We  are 
particularly  interested  in  aeroacoustic  tracking  of  ground  vehicles  using  an 
array  of  microphone  arrays.  Signal  coherence  is  known  to  degrade  with 
increased  sensor  separation  for  low  frequency  sounds  (10-300  Hz)  prop¬ 
agating  through  air  (for  examples,  see  Wilson  [7,  8]).  Thus  it  is  important 
to  understand  the  fundamental  limitations  of  TD  and  Doppler  estimation 
when  the  signals  are  -partially  coherent,  and  we  provide  a  detailed  study  of 
this  question  in  this  report.  Our  results  quantify  the  scenarios  in  which  TDE 
is  feasible  as  a  function  of  signal  coherence,  signal-to-noise  ratio  (SNR)  per 
sensor,  fractional  bandwidth  of  the  signal,  and  time-bandwidth  product  of 
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the  observed  data.  The  basic  result  is  that  for  a  given  SNR,  fractional  band¬ 
width,  and  time-bandwidth  product,  there  exists  a  "threshold  coherence" 
value  that  must  be  exceeded  in  order  for  TDE  to  achieve  the  Cramer-Rao 
Bound  (CRB).  The  analysis  is  based  on  Ziv-Zakai  bounds  for  TDE,  as  de¬ 
scribed  by  Weiss  and  Weinstein  [9,10]. 

The  specific  tasks  in  the  statement  of  work  for  this  contract  are  as  follows. 

a.  Study  alternative  approaches  to  modeling  source  motion  and  the  im¬ 
pact  on  localization  and  tracking  performance.  Two  possible 
approaches  are  as  follows: 

-  Model  the  sources  as  nonmoving  over  a  short  time  interval.  The 
validity  of  this  model  depends  on  the  source  velocity  and  the 
length  of  the  time  interval. 

-  Describe  the  source  motion  over  a  time  interval  with  a  paramet¬ 
ric  model.  Evaluate  the  tradeoffs  in  processing  complexity  and 
tracking  performance  associated  with  each  approach. 

b.  Investigate  the  performance  of  algorithms  for  TDE  with  partially  co¬ 
herent  signals  that  are  measured  at  widely  separated  aeroacoustic 
sensors. 

c.  Study  the  impact  of  compensating  for  the  Doppler  effect  that  is  caused 
by  the  source  motion.  Evaluate  the  performance  improvements 
and  processing  complexity  that  are  associated  with  Doppler 
compensation. 

d.  Investigate  performance  bounds  on  localization  accuracy  (Cramer- 
Rao,  Ziv-Zakai,  and/or  Barankin)  that  include  the  signal  structure 
and  the  partial  signal  coherence  at  distributed  arrays. 

e.  Investigate  algorithms  for  joint  tracking  and  identification  that  share 
information  about  source  location,  heading,  orientation,  and  source 
signal  attributes. 

This  task  was  elaborated  in  the  modification  of  July  2001  as  follows: 

Classification  of  ground  vehicles  is  currently  performed  using  fea¬ 
tures  derived  solely  from  the  signal  waveforms  measured  at  dis¬ 
tributed  sensors.  Classification  bounds  will  be  studied  that  quantify 
the  amount  of  additional  information  that  is  available  if  the  following 
data  are  used  in  conjunction  with  the  signal  waveforms:  target  aspect 
angle,  range,  and  heading /velocity  (obtained  from  the  tracking  algo¬ 
rithm),  and  coherence  level  of  the  signals  measured  at  widely  sepa¬ 
rated  sensors.  Acoustic  waveforms  for  ground  vehicles  as  a  function 
of  range,  aspect  angle,  and  heading /velocity  will  be  obtained  from 
existing  data  files  containing  moving  vehicles  and  distributed  sen¬ 
sors.  New  measurements  may  also  be  required  to  obtain  more  com¬ 
plete  data  sets. 
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f.  The  following  task  was  added  in  the  modification  of  July  2001: 

The  previous  work  has  focused  on  localization  of  a  single  source 
using  widely  separated  sensors  with  partial  spatial  coherence.  In¬ 
vestigate  extensions  of  the  single-source  analysis  and  techniques  to 
multiple-source  scenarios. 

Task  a  is  addressed  in  section  5  of  this  report.  Tasks  b-d  are  closely  related 
and  are  addressed  in  sections  2-5.  Task  f  is  addressed  in  section  6,  where 
we  present  an  example  with  two  moving  sources.  Task  e  relates  to  classi¬ 
fication  of  ground  vehicles  and  is  not  addressed  in  this  report.  We  have 
developed  a  theoretical  methodology  for  quantifying  the  amount  of  infor¬ 
mation  provided  by  additional  "features  for  classification  [11],  but  we  have 
not  applied  the  methodology  to  acoustical  signatures  of  ground  vehicles. 

This  report  is  organized  as  follows.  The  sensor  data  model  is  presented  in 
section  2  for  the  case  of  a  nonmoving  source.  TDE  with  partially  coherent 
signals  is  studied  in  section  3,  and  section  4  presents  a  subspace  algorithm 
for  narrowband  source  localization  with  distributed  arrays  and  partially 
coherent  signals.  The  sensor  data  model  is  extended  to  moving  sources  in 
section  5,  and  an  algorithm  is  proposed  for  tracking  with  an  array  of  arrays. 
Section  6  includes  an  example  based  on  measured  aeroacoustic  data  for  two 
moving  ground  vehicles.  The  report  concludes  in  section  7  with  a  summary 
and  concluding  remarks. 
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2.  Data  Model  for  a  Nonmoving  Source 


A  model  is  formulated  in  this  section  for  the  discrete-time  signals  received 
by  the  sensors  in  an  array  of  arrays.  To  begin,  suppose  a  single  nonmoving 
source  is  located  at  coordinates  {xs,ys)  in  the  {x,y)  plane,  and  consider 
H  arrays  that  are  distributed  in  the  same  plane,  as  illustrated  in  Figure  1. 
Each  array  h  e  {1, ...  ,H}  contains  Nh  sensors  and  has  a  reference  sensor 
located  at  coordinates  {xh,  yh)-  The  location  of  sensor  n  G  {1, . . . ,  Nh}  is  at 
{Xh  +  Nxhn,yh  +  ^yhn),  where  {Axhn,^yhn)  is  the  relative  location  with 
respect  to  the  reference  sensor.  If  c  is  the  speed  of  propagation,  then  the 
propagation  time  from  the  source  to  the  reference  sensor  on  array  h  is 

rh  =  —  =  -  [(2:5  -  Xh)'^  +  {ys  -  yhf]  ,  (1) 

c  c 

where  dh  is  the  distance  from  the  source  to  array  h.  We  model  the  wave- 
fronts  over  individual  array  apertures  as  perfectly  coherent  plane  waves. 
Then  in  the  far-field  approximation,  the  propagation  time  from  the  source 
to  sensor  n  on  array  h  is  expressed  by  Th  +  Thn,  where 

1 

'Thn  ~ 

C 

is  the  propagation  time  from  the  reference  sensor  on  array  h  to  sensor  n  on 
array  h,  and  (/>/,  is  the  bearing  of  the  source  with  respect  to  array  h.  Note  that 
while  the  far-field  approximation  (2)  is  reasonable  over  individual  array 
apertures,  the  wavefront  curvature  that  is  inherent  in  (1)  must  be  retained 
in  order  to  model  wide  separations  between  arrays. 

The  time  signal  received  at  sensor  n  on  array  h  due  to  the  source  will 
be  represented  as  Sh{t  —  Th  —  Thn),  where  the  vector  of  signals  s{t)  = 
[si{t), . . .  received  at  the  H  arrays  are  modeled  as  real-valued, 

continuous-time,  zero-mean,  jointly  wide-sense  stationary,  Gaussian  ran¬ 
dom  processes  with  —00  <  t  <  00.  These  processes  are  fully  specified  by 
the  H  X  H  cross-correlation  function  matrix 

Rs{t)  =  E{s{t  +  T)s{t)'^},  (3) 


xs-xh  .  .  ys-yh  . 

AXhn  T  ,  Ayhn 


dh 


dh 


=  -  -  [{cos  (l>h)  Ax hn  +  {sin  (l)h)  Ay hn]  (2) 


Figure  1.  Geometry  of 
nonmoving  source 
location  and  an  array  of 
arrays.  A  communi¬ 
cation  link  is  available 
between  each  array  and 
the  fusion  center. 
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where  E  denotes  expectation,  superscript  T  denotes  transpose,  and  we  will 
later  use  the  notation  superscript  *  and  superscript  H  to  denote  complex 
conjugate  and  conjugate  transpose,  respectively.  The  {g,  h)  element  in  (3)  is 
the  cross-correlation  function 


rs,gh{T)  =  E{sg{t  -h  r)  Sh{t)}  (4) 

between  the  signals  received  at  arrays  g  and  h.  The  correlation  functions  (3) 
and  (4)  are  equivalently  characterized  by  their  Fourier  transforms,  which 
are  the  cross-spectral  density  functions 

/OO 

rs,gh{T)  exp(-jwT)  dr  (5) 

-OO 

and  the  associated  cross-spectral  density  matrix 

Gs{cv)  =  E{Ks{r)}.  (6) 

The  diagonal  elements  Gg^hhi^^)  of  (6)  are  the  power  spectral  density  (PSD) 
functions  of  the  signals  Sh{t),  and  hence  they  describe  the  distribution  of 
average  signal  power  with  frequency  The  model  allows  the  PSD  to  vary 
from  one  array  to  another  to  reflect  propagation  differences  and  source  as¬ 
pect  angle  differences. 

The  off-diagonal  elements  of  (6),  Gg^ghioo),  are  the  cross-spectral  density 
(CSD)  functions  for  the  signals  Sg{t)  and  Sh{t)  received  at  distinct  arrays 
g  ^  h.  In  general,  the  CSD  functions  have  the  form 

Gs,gh{‘^)  =  'ls,gh{‘^)  [Gs,gg{uj)Gs,hh{oj)]^^'^  ,  (7) 

where  ^s,gh{^)  is  the  spectral  coherence  function  for  the  signals,  which  has 
the  property  0  <  \^s,gh{co)\  <  1.  Coherence  magnitude  \  js,gh{‘^)\  =  1  corre¬ 
sponds  to  perfect  correlation  between  the  signals  at  sensors  g  and  h,  while 
the  partially  coherent  case  \^s,gh{^)  \  <  1  models  random  effects  in  the  prop¬ 
agation  paths  from  the  source  to  sensors  g  and  h.  Note  that  our  assumption 
of  perfect  spatial  coherence  across  individual  arrays  implies  that  the  ran¬ 
dom  propagation  effects  have  negligible  impact  on  the  intra-array  delays 
Thn  in  (2)  and  the  bearings  0i, . . .  fn- 

The  signal  received  at  sensor  n  on  array  h  is  modelled  as  a  sum  of  the 
delayed  source  signal  and  noise, 

Zhn{t)  =  Sh{t  -Th-  Thn)  +  VJhn{t),  (8) 


where  the  noise  signals  Whn{t)  are  modeled  as  real-valued,  continuous¬ 
time,  zero-mean,  jointly  wide-sense  stationary,  Gaussian  random  processes 
that  are  mutually  uncorrelated  at  distinct  sensors,  and  are  uncorrelated 
from  the  signals.  That  is,  the  noise  correlation  properties  are 

E{Wgm{t  -h  T)whn{t)}  =  r^(r)  Sghdmn  (9) 


and 

E{Wgm{t  +  T)Sh{t)}  =  0,  (10) 
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where  r^(r)  is  the  noise  autocorrelation  function,  and  the  noise  PSD  is 
Gw{u!)  =  JP{ru,(r)}.  We  then  collect  the  observations  at  each  array  h  into 
NhX  1  vectors  Zh{t)  =  [zhi{t), Zh,Nh{'t)]'^  for  h  =  1, . . . ,  H,  and  we  fur¬ 
ther  collect  the  observations  from  the  H  arrays  into  a  (A^i  -I-  •  •  •  -I-  Nh)  x  1 
vector 

Zi(f) 

m  =  ;  •  (11) 

The  elements  of  Z{t)  in  (11)  are  zero-mean,  jointly  wide-sense  stationary, 
Gaussian  random  processes.  We  can  express  the  CSD  matrix  of  Z(f)  in  a 
convenient  form  with  the  following  definitions.  The  array  manifold  for  ar¬ 
ray  h  at  frequency  uj  is 

exp{-juThi)  exp  {{cos  (ph)  Ax  hi  +  {sin  (jih)  Ay  hi)] 

nh{u)  =  :  =  :  ,  (12) 

_  exp{-jujTh,Nh)  \  I  exp  [j^  ((cos  (f)h)Axh,Nh  +  (sin  (i)h)Ayh,Nh)]  . 

using  Thn  from  (2)  and  assuming  that  the  sensors  have  omnidirectional  re¬ 
sponse  to  sources  in  the  plane  of  the  array.  Let  us  define  the  relative  TD  of 
the  signal  at  arrays  g  and  h  as 

Dgh  =  rg-  Th,  (13) 

where  r/,  is  defined  in  (1).  Then  the  cross-spectral  density  matrix  of  Z{t)  in 
(11)  has  the  form 

Gz(ca)  = 

(14) 

ai{uj)ai{uj)^Gs,ii{uj)  ■■■  ai{u;)aH{uj)^  exp{-jujDiH)Gs,iH{uj) 

;  '  •  _  :  -|-  Gw{u!)1. 

aH{uj)ai{uj)^  exp{+j(jjDiH)Gs,iH{>^)*  •••  aH{u;)aH{uj)^  Gs,hh{uj) 

The  source  CSD  functions  Gs,gh{^)  in  (14)  can  be  expressed  in  terms  of  the 
signal  spectral  coherence  ^s,gh{'^)  using  (7).  Note  that  (14)  depends  on  the 
source  location  parameters  {xs,ys)  through  the  bearings  (j)h  in  a/,(w)  and 
the  pairwise  TD  differences  Dgh- 

2.1  Cramer-Rao  Bound  (CRB) 

The  CRB  provides  a  lower  bound  on  the  variance  of  any  unbiased  es¬ 
timator.  The  problem  of  interest  is  estimation  of  the  source  location  pa¬ 
rameter  vector  0  =  [xs,ys]'^  using  T  samples  of  the  sensor  signals 
Z(0),  Z(Ts), . . . ,  Z((T  —  1)  •  Ts),  where  Tg  is  the  sampling  period.  The  total 
observation  time  is  T  =  T-Ts.  Let  us  denote  the  sampling  rate  by  =  l/Tg 
and  iOs  =  “^nfs.  We  will  assume  that  the  continuous-time  random  processes 
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Z(t)  are  band-limited,  and  that  the  sampling  rate  fs  is  greater  than  twice 
the  bandwidth  of  the  processes.  Then  Friedlander  and  Whittle  [12, 13]  have 
shown  that  the  Fisher  information  matrix  (FIM)  J  for  the  parameters  0 
based  on  the  samples  Z(0),  Z{Ts), . . . ,  Z((r  —  1)  •  Tg)  has  elements 


j.-L 

~  dvr 


tri^^Gz 

d6i 


(  \  — ^  /  \  —  l 


duj,  i,j  =  1,2, 


(15) 


where  "tr"  denotes  the  trace  of  the  matrix.  The  CRB  matrix  C  =  then 
has  the  property  that  the  covariance  matrix  of  any  unbiased  estimator  0 
satisfies  Cov(0)  —  C  >  0,  where  >  0  means  that  Cov(0)  —  C  is  positive 
semidefinite  [14].  Equation  (15)  provides  a  convenient  way  to  compute  the 
FIM  for  the  array  of  arrays  model  as  a  function  of  the  signal  coherence  be¬ 
tween  distributed  arrays,  the  signal  and  noise  bandwidth  and  power  spec¬ 
tra,  and  the  sensor  placement  geometry.  The  CRB  is  evaluated  for  various 
scenarios  in  our  previous  work  [6],  and  additional  examples  are  presented 
in  section  2.2.1  of  this  report. 

Let  us  consider  the  CRB  for  an  acoustic  source  that  has  a  narrowband  power 
spectrum,  i.e.,  the  PSD  Gg^hhi^)  of  the  signal  at  each  array  h  =  1, ...  ,H 
is  nonzero  only  in  a  narrow  band  of  frequencies  loq  —  (Aca/2)  <  u)  < 
wq  +  {A(jj/2).  If  the  bandwidth  Aw  is  chosen  small  enough  so  that  the  w- 
dependent  quantities  in  (15)  are  well  approximated  by  their  value  at  loq, 
then  the  narrowband  approximation  to  the  FIM  (15)  is 


J, 


TAlo 


-tr 


de^ 


de 


3 


(16) 


The  quantity  multiplying  the  FIM  in  (16)  is  the  time-bandwidth  prod¬ 
uct  of  the  observations.  In  narrowband  array  processing,  the  T  time  sam¬ 
ples  per  sensor  are  often  segmented  into  M  blocks  containing  T/M  sample 
each.  Then  the  discrete  Fourier  transform  (DFT)  is  applied  to  each  block, 
and  the  complex  coefficients  at  frequency  wq  (at  each  sensor)  are  used  to 
form  M  array  "snapshots."  In  this  case,  the  quantity  is  approximately 
equal  to  M,  the  number  of  snapshots. 


2.2  Examples 

Next,  we  present  numerical  examples  that  evaluate  the  CRB  on  localization 
accuracy  in  (15)  for  a  narrowband  and  a  wideband  source.  We  also  show 
measured  aeroacoustic  data  from  a  ground  vehicle  and  evaluate  the  spec¬ 
tral  coherence  of  the  source  at  widely  separated  sensors. 

2.2.1  CRB  Evaluation 

In  our  previous  work  [6],  we  developed  CRBs  on  localization  accuracy  for 
the  following  three  schemes. 

1.  Each  array  estimates  the  source  bearing  and  transmits  the  bearing  es¬ 
timate  to  the  fusion  center.  The  fusion  processor  then  triangulates  the 
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bearings  to  estimate  the  source  location.  This  approach  does  not  ex¬ 
ploit  wavefront  coherence  between  the  distributed  arrays,  but  it  does 
minimize  the  communication  bandwidth  required  to  transmit  data 
from  the  arrays  to  the  fusion  center. 

2.  The  raw  data  from  all  sensors  is  jointly  processed  to  estimate  the 
source  location.  This  is  the  optimum  approach  and  it  fully  utilizes  the 
coherence  between  distributed  arrays.  However,  it  requires  a  large 
communication  bandwidth,  since  the  data  from  all  of  the  sensors  must 
be  transmitted  to  the  fusion  center. 

3.  Combination  of  methods  1  and  2:  The  objective  is  to  perform  some 
processing  at  the  individual  arrays  to  reduce  the  communication  band¬ 
width  requirement  while  still  exploiting  the  coherence  between  dis¬ 
tributed  arrays.  Each  array  estimates  the  source  bearing  and  trans¬ 
mits  the  bearing  estimate  to  the  fusion  center.  In  addition,  the  raw 
data  from  one  sensor  in  each  array  is  transmitted  to  the  fusion  cen¬ 
ter.  The  fusion  center  then  estimates  the  propagation  time  delay  be¬ 
tween  pairs  of  distributed  arrays,  and  triangulates  these  time  delay 
estimates  with  the  bearing  estimates  to  localize  the  source. 

The  first  scheme  is  "incoherent  triangulation"  of  the  bearing  estimates,  and 
it  is  obtained  from  the  CRB  in  (15)  by  setting  the  coherence  parameters  to 
zero,  7g/i  =  0,  in  (7)  and  (14).  The  second  scheme  is  obtained  from  the  CRB 
in  (15)  with  nonzero  coherence,  and  the  CRB  for  the  third  scheme  is 
described  in  [6].  Our  objective  in  this  section  is  to  evaluate  these  CRBs  for 
a  narrowband  source  and  a  wideband  source.  We  will  refer  back  to  these 
examples  in  section  3  after  developing  the  threshold  coherence  analysis  for 
TDE.  The  analysis  in  section  3  will  show  that  the  CRBs  for  the  narrowband 
source  case  are  unachievable  when  there  is  any  loss  of  signal  coherence  be¬ 
tween  arrays,  while  the  CRBs  for  the  wideband  source  case  are  achievable 
when  moderate  coherence  losses  occur. 

Consider  H  =  3  identical  arrays,  each  of  which  contains  Ni  =  ■  ■  ■  =  Nh  = 

7  sensors.  Each  array  is  circular  with  4-ft  radius,  with  six  sensors  equally 
spaced  around  the  perimeter  and  one  sensor  in  the  center.  We  first  evalu¬ 
ate  the  CRB  for  a  narrowband  source  with  a  1-Hz  bandwidth  centered  at 
50  Hz  and  SNR  =  10  dB  at  each  sensor.  That  is,  Gs,hh{^)/Gw{u!)  =  10  for 
h  =  and  27r(49.5)  <  u)  <  27r(50.5)  rad/s.  The  signal  coherence 

ls,gh{^)  =  7^(1^)  is  varied  between  0  and  1.  We  assume  that  T  =  4000  time 
samples  are  obtained  at  each  sensor  with  sampling  rate  /«  =  2000  sam¬ 
ples/  s.  The  source  localization  performance  is  evaluated  by  computing  the 
ellipse  in  (x,  y)  coordinates  that  satisfies  the  expression 


X 

X  y 

J 

.  y . 

(17) 


where  J  is  the  EIM  in  (15).  If  the  errors  in  (x,  y)  localization  are  jointly  Gaus¬ 
sian  distributed,  then  the  ellipse  (17)  represents  the  contour  at  one  stan¬ 
dard  deviation  in  root-mean-square  (RMS)  error.  The  error  ellipse  for  any 
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unbiased  estimator  of  source  location  cannot  be  smaller  than  this  ellipse 
derived  from  the  FIM. 

The  H  =  3  arrays  are  located  at  coordinates  {xi,yi)  =  (0,0),  (x2,y2)  = 
(400, 400),  and  (xs,  ys)  =  (100, 0).  One  source  located  at  (xs,  y*)  =  (200, 300), 
where  the  units  are  meters,  as  illustrated  in  Figure  2a.  The  RMS  error  el¬ 
lipses  for  joint  processing  of  all  sensor  data  for  coherence  values  7s (w)  = 
0, 0.5,  and  1  are  also  shown  in  Figure  2a.  The  coherence  between  all  pairs 
of  arrays  is  assumed  to  be  identical,  i.e.,  7s,g/i(w)  =  7s(tc)  for  {g,h)  = 
(1,2),  (1,3),  (2, 3).  The  largest  ellipse  in  Figure  2a  corresponds  to  incoherent 
signals  7s  (ta)  =  0  and  characterizes  the  performance  of  the  simple  method 
of  triangulation  using  the  bearing  estimates  from  the  three  arrays.  Figure  2b 
shows  the  ellipse  radius  for  various  values  of  the  signal  coherence  7s  (a;) . 
The  ellipses  for  75(0;)  =  0.5  and  1  are  difficult  to  see  in  Figure  2a  because 
they  fall  on  the  lines  of  the  x  that  marks  the  source  location,  illustrating 
that  signal  coherence  between  the  arrays  significantly  improves  the  CRB 
on  source  localization  accuracy.  Note  also  that  for  this  scenario,  the  local¬ 
ization  scheme  based  on  bearing  estimation  with  each  array  and  TDE  using 
one  sensor  from  each  array  performs  equivalently  to  the  optimum  joint  pro¬ 
cessing  scheme.  Figure  2c  shows  a  closer  view  of  the  error  ellipses  for  the 
scheme  of  bearing  estimation  plus  TDE  with  one  sensor  from  each  array. 
The  ellipses  are  identical  to  those  in  Figure  2a  for  joint  processing. 


Figure  2.  RMS  source 
localization  error 
ellipses  based  on  the 
CRB  for  H  =  3  arrays 
and  one  narrowband 
source.  The  array  and 
source  locations  are 
shown  in  (a),  along  with 
the  error  ellipses  for 
joint  processing  of  all 
sensor  data  for 
coherence  values 


CRB  Ellispes  for  Coherence  0,  0.5, 1 .0  (Joint  Processing) 
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7s(a;)  =  0, 0.5,  and  1. 
Part  (b)  shows  the  error 
ellipse  radius  (major 
axis^  +  minor  axis^)^/^ 
for  a  range  of  coherence 
values,  comparing  joint 
processing  with  the 
reduced-complexity 
scheme  of  bearing 
estimation  plus  TDE 
using  data  from  one 
sensor  per  array.  Part 
(c)  is  a  closer  view  of  the 
RMS  error  ellipses  for 
the  bearing  plus  TDE 
scheme. 
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Figure  3  presents  corresponding  results  for  a  wideband  source  with  band¬ 
width  20  Hz  centered  at  50  Hz  and  SNR  16  dB.  That  is,  Gs,hh/Gw  =  40 
for  27r(40)  <  to  <  27r(60)  rad/s,  h  =  1,. . .  ,H.  T  =  2000  time  samples  are 
obtained  at  each  sensor  with  sampling  rate  /*  =  2000  samples/s,  so  the 
observation  time  is  1  s.  As  in  the  narrowband  case  in  Figure  2,  Figure  3 
shows  that  joint  processing  reduces  the  CRB  compared  with  bearings-only 
triangulation,  and  bearing  plus  TD  estimation  is  nearly  optimum. 

The  CRB  provides  a  lower  bound  on  the  variance  of  unbiased  estimates, 
so  an  important  question  is  whether  an  estimator  can  achieve  the  CRB.  We 
show  in  section  3  that  the  coherent  processing  CRBs  for  the  narrowband 
scenario  illustrated  in  Figure  2  are  achievable  only  when  the  coherence  is 
perfect,  i.e.,  =  1.  Therefore  for  that  scenario,  bearings-only  triangulation 

is  optimum  in  the  presence  of  even  small  coherence  losses.  However,  for  the 
wideband  scenario  illustrated  in  Figure  2,  the  coherent  processing  CRBs  are 
achievable  for  coherence  values  >  0.75,  so  the  wideband  source  provides 
tolerance  to  moderate  values  of  coherence  loss. 

2.2.2  Coherence  in  Measured  Data 

Next,  we  present  results  from  measured  aeroacoustic  data  to  illustrate  typi¬ 
cal  values  of  signal  coherence  at  distributed  arrays.  The  experimental  setup 
is  illustrated  in  Figure  4a,  which  shows  the  path  of  a  moving  ground  vehi¬ 
cle  and  the  locations  of  four  microphone  arrays  (labeled  1,  3,  4,  5).  Each 
array  is  circular  with  N  =  7  sensors,  4-ft  radius,  and  six  sensors  equally 


Figure  3.  RMS  source 
localization  error 
ellipses  based  on  the 
CRB  for  H  —  S  arrays 
and  one  wideband 
source.  The  array  and 
source  locations  are 
shown  in  (a),  along  with 
the  error  ellipses  for 
joint  processing  of  all 
sensor  data  for 
coherence  values 
ys(w)  =  0,  0.5,  and  1. 
Part  (b)  shows  the  error 
ellipse  radius, 
comparing  joint 
processing  with  the 
reduced-complexity 
scheme  of  bearing 
estimation  plus  TDE 
using  data  from  one 
sensor  per  array.  Part 
(c)  is  a  closer  view  of  the 
RMS  error  ellipses  for 
the  bearing  plus  TDE 
scheme. 
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spaced  around  the  perimeter  with  one  sensor  in  the  center.  We  focus  on  the 
10-s  segment  indicated  by  the  in  Figure  4a  (which  correspond  to  the 
time  segment  340-350  s  in  the  data).  Figure  4b  shows  the  power  spectral 
density  (PSD)  of  the  data  measured  at  arrays  1  and  3  during  the  10-s  seg¬ 
ment.  Note  the  dominant  harmonic  at  39  Hz.  Figure  4c  shows  the  estimated 
coherence  between  arrays  1  and  3  during  the  10-s  segment.  The  coherence 
is  approximately  0.85  at  40  Hz,  which  demonstrates  the  presence  of  sig¬ 
nificant  coherence  at  widely  separated  microphones.  Figure  4c  shows  the 
estimated  coherence  between  two  sensors  on  array  1,  spaced  by  8  ft.  Note 
that  the  coherence  is  close  to  unity  for  frequencies  in  the  range  from  about 
40  to  200  Hz,  so  our  model  of  perfect  signal  coherence  over  individual  ar¬ 
rays  seems  reasonable. 


Figure  4.  (a)  Path  of 
ground  vehicle  and 
array  locations  for 
measured  data. 

(b)  Mean  PSD  at  arrays 
1  and  3  estimated  from 
measured  data  over  the 
10-s  segment  0  in  (b). 
Top  panel  is  Gs,  ii(/)/ 
bottom  panel  is 
G«,33(/).(c)Mean 
spectral  coherence 
7s,i3(/)  between  arrays 
1  and  3  estimated  over 
the  10-s  segment. 

(d)  Mean  spectral 
coherence  for  two 
sensors  on  array  1,  with 
sensor  spacing  8  ft. 


(a)  Ground  Vehicle  Path  and  Array  Locations 
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The  Doppler  effect  due  to  source  motion  was  compensated  prior  to  the  co¬ 
herence  estimate  shown  in  Figure  4c.  Without  Doppler  compensation,  the 
coherence  is  significantly  reduced,  as  shown  in  Figure  5a.  The  time-varying 
radial  velocity  of  the  source  with  respect  to  each  array  in  Figure  4a  is  plot¬ 
ted  in  the  top  panel  of  Figure  5b.  If  s{t)  is  the  waveform  emitted  by  the 
source  that  is  moving  with  radial  velocity  v  with  respect  to  the  sensor,  then 
the  sensor  receives  a  waveform  with  the  form  s{at),  where  the  scaling  fac¬ 
tor  a  is 

0  =  1--  (18) 
c 

and  c  is  the  speed  of  wave  propagation.  The  scaling  factor  a  is  plotted  in 
the  bottom  panel  of  Figure  5b.  Note  that  for  this  data  set,  0.98  <  a  <  1.02, 
which  corresponds  to  a  Doppler  frequency  shift  of  approximately  ±1  Hz 
for  an  emitted  tone  at  50  Hz.  We  use  a  digital  resampling  algorithm  to  com¬ 
pensate  for  the  Doppler  effect. 


Figure  5.  (a)  Mean 
short-time  spectral 
coherence  between 
arrays  1  and  3  if 
Doppler  is  not 
compensated,  (b)  Radial 
velocity  and  Doppler 
scaling  factor  a  in  (18) 
for  source  in  Figure  4a 
with  respect  to  each 
array. 
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3.  Time  Delay  Estimation  (TDE) 


CSD 
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e 


Figure  6.  TDE  problem 
for  a  nonmoving  source 
with  H  =  2  sensors. 


The  CRB  results  presented  in  section  2.2.1  indicate  that  TDE  between  widely 
spaced  sensors  is  an  effective  way  to  improve  the  source  localization  ac¬ 
curacy  with  joint  processing.  Fundamental  performance  limits  for  passive 
time  delay  and  Doppler  estimation  have  been  studied  extensively  for  sev¬ 
eral  decades  (e.g.,  see  the  collection  of  papers  edited  by  Carter  [15]).  The 
fundamental  limits  are  usually  parameterized  in  terms  of  the  SNR  at  each 
sensor,  the  spectral  support  of  the  signals  (fractional  bandwidth),  and  the 
time-bandwidth  product  of  the  observations.  We  are  particularly  interested 
in  the  application  of  aeroacoustic  tracking  of  ground  vehicles  using  an  ar¬ 
ray  of  microphone  arrays.  The  coherence  between  the  signals  is  known  to 
degrade  with  increased  spatial  separation  between  the  signals  [7,  8]  due  to 
random  motion  of  the  air.  This  coherence  loss  significantly  affects  the  time 
delay  and  Doppler  estimation  accuracy. 

In  order  to  quantify  the  effect  of  partial  signal  coherence  on  TD  and 
Doppler  estimation,  we  present  Cramer-Rao  and  Ziv-Zakai  bounds  that  are 
explicitly  parameterized  by  the  signal  coherence,  along  with  the  traditional 
parameters  of  SNR,  fractional  bandwidth,  and  time-bandwidth  product. 
The  basic  result  is  that  for  a  given  SNR,  fractional  bandwidth,  and  time- 
bandwidth  product,  there  exists  a  "threshold  coherence"  value  that  must 
be  exceeded  in  order  for  TDE  to  achieve  the  CRB.  The  analysis  is  based  on 
Ziv-Zakai  bounds  for  TDE,  as  in  Weiss  and  Weinstein  [9, 10]. 

Let  us  parameterize  the  model  in  (14)  by  the  bearings  cph  and  the  time-delay 
differences  Dgh,  and  consider  first  the  case  of  H  =  2  sensors.  Figure  6  con¬ 
tains  an  illustration  of  the  TDE  problem,  with  the  differential  TD  defined 
as  D  =  D21.  It  follows  from  (14)  that  the  CSD  matrix  of  the  sensor  data  in 
Figure  6  is 


GzH  = 


(19) 


1^8,11(1^)  +  Gyj{U}) 
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The  signal  coherence  function  7^,12 (<^)  describes  the  degree  of  correlation 
that  remains  in  the  signal  emitted  by  the  source  at  each  frequency  uj  af¬ 
ter  propagating  to  sensors  1  and  2.  Next,  we  develop  an  SNR-like  expres¬ 
sion  for  the  two-sensor  case  that  appears  in  all  subsequent  expressions 
for  fundamental  limits  on  TD  and  Doppler  estimation.  We  begin  with  the 
magnitude-squared  coherence  (MSC)  [15]  of  the  observed  signals  zi{t),Z2{t) 
as  a  function  of  the  signal  coherence  magnitude,  |7s,i2(<^)|/  and  other  spec¬ 
tral  density  parameters: 


MSC,(|7.,i2MI) 


|CSD[2;i(f),2;2(f)]P  _  Gs,Ii{uj)Gs,22{^) 

PSD[zi(f)]  •  PSD[z2(t)]  “  [G,,ii(u;)  +  GM]  [GsM^)  + 


l7s,i2(w)P 


Y  Gw 


-1 


1  + 


/  Gs,22(‘^) 
Y  Gw{^) 


-1 


<  1 . 


(20) 


Then  the  following  "SNR"  expression  appears  in  subsequent  performance 
bounds: 


SNR(|7,,i2(a;)|) 


MSC,(|7«,i2(u;)|) 

1-MSC,(|7m2MI) 


< 


\|75,12(w)I^  |_ 
l7s,12(^)P 
1  -  |7s,12(‘^)l^ 


(Gs,n{^) 
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(21) 

l|  (22) 
(23) 


The  inequality  (23)  shows  that  signal  coherence  loss  (|7s,i2(<^)  |  <  1)  severely 
limits  the  "SNR"  quantity  that  characterizes  performance,  even  if  the  SNR 
per  sensor  Gs,ii{oj)/Gw{u!)  is  very  large. 


3.1  Bounds  for  TDE 


We  can  use  (19)  in  (15)  to  find  the  CRB  for  TDE  with  H  =  2  sensors,  yielding 


CRB(D) 


47r 

Y 


-1 


ta^SNR(|7,j^i2(t(;)|)  doj 


(24) 


where  T  is  the  total  observation  time  of  the  sensor  data  and  SNR  (|7s,i2(:^)|) 
is  defined  in  (22).  Let  us  consider  the  case  in  which  the  signal  PSDs,  the 
noise  PSD,  and  the  coherence  are  flat  (constant)  over  a  bandwidth  Aca  rad/ s 
centered  at  ujq  rad/ s.  If  we  omit  the  frequency  dependence  of  0^,22,  G^, 

and  7s, 12/  then  the  integral  in  (24)  may  be  evaluated  to  yield  the  following 
CRB  expression: 
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The  quantity  is  the  time-bandwidth  product  of  the  observations, 

(^)  is  the  fractional  bandwidth  of  the  signal,  and  Gg^hh/Gw  is  the  SNR 
at  sensor  h.  Note  from  the  high-SNR  limit  in  (26)  that  when  the  signals  are 
partially  coherent  1 75,12!  <  1/  increased  source  power  does  not  reduce  the 
CRB.  Improved  TDE  accuracy  is  obtained  with  partially  coherent  signals  by 
increasing  the  observation  time  T  or  changing  the  spectral  support  of  the 
signal,  which  is  [wq  —  /2 ,  wq  -|-  Aca /2] .  The  spectral  support  of  the  signal  is 
not  controllable  in  passive  TDE  applications,  so  increased  observation  time 
is  the  only  means  for  improving  the  TDE  accuracy  with  partially  coherent 
signals.  Source  motion  becomes  more  important  during  long  observation 
times,  and  in  section  5  we  extend  the  model  to  include  source  motion. 

The  CRB  in  (25)  agrees  with  known  results  for  perfectly  coherent  signals 
[15]  and  with  results  from  the  medical  ultrasound  literature  [16,  17]  for 
partially  correlated  speckle  signals.  The  medical  ultrasound  application  is 
distinguished  from  the  aeroacoustic  tracking  of  ground  vehicles  in  that  the 
former  is  typically  an  active  system  while  the  latter  is  passive.  The  medical 
ultrasound  application  therefore  allows  much  more  control  over  the  SNR 
and  bandwidth  of  the  signals.  In  passive  aeroacoustics,  the  received  signals 
are  emitted  by  a  vehicle  and  are  not  controllable  for  the  purposes  of  TD 
estimation. 

With  perfectly  coherent  signals,  it  is  well  known  that  the  CRB  on  TDE  is 
achievable  only  when  the  SNR  expression  in  (22)  (with  175,12(1^)1  =  1)  ex¬ 
ceeds  a  threshold  [9, 10].  Next  we  show  that  for  TDE  with  partially  coherent 
signals,  a  similar  threshold  phenomenon  occurs  with  respect  to  coherence. 
That  is,  the  coherence  must  exceed  a  threshold  in  order  to  achieve  the  CRB 
(24)  on  TDE.  We  state  the  threshold  coherence  formula  for  the  following 
simplified  scenario.  The  signal  and  noise  spectra  are  flat  over  a  bandwidth 
of  Au!  rad/s  centered  at  tco  rad/s,  and  the  observation  time  is  T  seconds. 
Eurther,  assume  that  the  signal  PSDs  are  identical  at  each  sensor,  and  define 
the  following  constants  for  notational  simplicity: 

^^,11(^0)  =  6^5, 22(1^0)  =  Gg,  Gw{iVo)  =  Gw,  and  75,12(^0)  =  Is-  (27) 
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Then  the  SNR  expression  in  (22)  has  the  form 


SNR(|7.|) 


1 


{Gs/G^ 


-1 


-  1 


(28) 


The  Ziv-Zakai  bound  developed  by  Weiss  and  Weinstein  [9, 10]  shows  that 
the  threshold  SNR  for  CRB  attainability  is  a  function  of  the  time-bandwidth 
product  {^2^)  and  the  fractional  bandwidth 


SNRthresh 


.-1 


(29) 


where  (p(y)  =  l/\/^  exp(— 1^/2)  dt.  It  follows  that  the  threshold  coher¬ 

ence  value  is 
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For  a  specific  TDE  scenario,  the  threshold  SNR  for  CRB  attainability  is 
given  by  (29),  and  (30)  provides  a  corresponding  threshold  coherence.  Since 
|7sP  <  1,  (30)  is  useful  only  if  Gs/G^  >  SNRthresh- 

Figure  7  contains  a  plot  of  (30)  for  a  particular  case  in  which  the  signals  are 
in  a  band  centered  at  ujq  =  27r50  rad/s  and  the  time  duration  is  T  =  2  s. 
Figure  7  shows  the  variation  in  threshold  coherence  as  a  function  of  sig¬ 
nal  bandwidth,  Ata.  Note  that  nearly  perfect  coherence  is  required  when 
the  signal  bandwidth  is  less  than  5  Hz  (or  10%  fractional  bandwidth).  The 
threshold  coherence  drops  sharply  for  values  of  signal  bandwidth  greater 
than  10  Hz  (20%  fractional  bandwidth).  Thus  for  sufficiently  wideband  sig¬ 
nals,  e.g..  Aw  >  27rl0  rad/s,  a  certain  amount  of  coherence  loss  can  be 
tolerated  while  still  allowing  unambiguous  TDE. 

Figures  8a-c  contain  plots  of  the  threshold  coherence  in  (30)  as  a  function 
of  the  time-bandwidth  product  {^2^),  SNR  and  fractional  bandwidth 

(^).  Note  that  ^  =  10  dB  is  nearly  equivalent  to  ^  ^  00.  We  note 
that  very  large  time-bandwidth  product  is  required  to  overcome  coherence 
loss  when  the  fractional  bandwidth  is  small  at  0.1.  The  variation  of  thresh¬ 
old  coherence  with  fractional  bandwidth  is  illustrated  in  Figure  8d.  For  a 
fixed  threshold  coherence  value,  such  as  0.7,  each  doubling  of  the  fractional 
bandwidth  reduces  the  required  time-bandwidth  product  by  about  a  factor 
of  10. 

Let  us  examine  a  narrowband  signal  scenario  that  is  typical  in  aeroacous- 
tics,  with  center  frequency  /o  =  a;o/(27r)  =  50  Hz  and  bandwidth  A/  = 
Aw/(27r)  =  5  Hz,  so  the  fractional  bandwidth  is  A///o  =  0.1.  From  Fig¬ 
ure  8a,  signal  coherence  =  0.8  requires  time-bandwidth  product  A/  • 
T  >  200,  so  the  necessary  time  duration  T  =  40  s  for  TDE  may  be  imprac¬ 
tical  for  moving  sources. 
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Figure  7.  Threshold 
coherence  vs. 
bandwidth  based  on 
(30)  for  loq  =  27r50  rad/ s 
and  GsjGw  =  0, 10,  and 
oo  dB. 


Figure  8.  Threshold 
coherence  value  from 
(30)  vs.  fime-bandwidfh 
producf  and 

SNR  G./Gjior 
fracfional  bandwidfh 
'values  (^)  (a)  0.1, 

(b)  0.5,  and  (c)  1.0.  In 
(d),  fhe  high  SNR 
curves  Gg/Gw  oo  are 
superimposed  for 
several  values  of 
fracfional  bandwidfh. 


(Og  =  2  71  50  rad/s,  T  =  2  s 


Larger  time-bandwidth  products  of  the  observed  signals  are  required  in 
order  to  make  TDE  feasible  in  environments  with  signal  coherence  loss.  As 
discussed  previously,  only  the  observation  time  is  controllable  in  passive 
applications,  thus  leading  us  to  consider  source  motion  models  in  section  5 
for  use  during  long  observation  intervals. 
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We  can  evaluate  the  threshold  coherence  for  the  narrowband  and  wideband 
scenarios  considered  in  section  2.2.1  for  the  CRB  examples.  The  results  are 
as  follows,  using  (29)  and  (30): 

•  Narrowband  case:  Gs/Gw  =  10,  ujq  =  27r50  rad/s,  Aca  =  27r  rad/s, 
T  =  2s 

Threshold  coherence  1. 

•  Wideband  case:  Gg/Gw  =  40,  tco  =  27r50  rad/s.  Ate  =  27r  •  20  rad/s, 
T  =  Is 

Threshold  coherence  0.75. 

Therefore,  for  the  narrowband  case,  joint  processing  of  the  data  from  dif¬ 
ferent  arrays  will  not  achieve  the  CRBs  in  Figure  2  when  there  is  any  loss  in 
signal  coherence.  For  the  wideband  case,  joint  processing  can  achieve  the 
CRBs  in  Figure  3  for  coherence  values  that  exceed  0.75. 

The  remainder  of  this  section  continues  to  focus  on  nonmoving  sources, 
with  a  simulation  example  presented  in  section  3.2  that  verifies  the  CRB 
and  threshold  coherence  values  for  TDE.  In  section  3.3,  the  H  =  2  sensor 
case  of  this  subsection  is  extended  to  TDE  with  H  >  2  sensors.  Section  3.4 
contains  examples  of  TDE  with  measured  aeroacoustic  data. 

3.2  TDE  Simulation  Examples 

Consider  TDE  at  H  =  2  sensors  with  varying  signal  coherence  Our  first 
simulation  example  involves  a  signal  with  reasonably  wide  bandwidth, 
A/  =  30  Hz  centered  at  /o  =  100  Hz,  so  the  fractional  bandwidth  A///o  = 
0.3.  The  signal,  noise,  and  coherence  are  flat  over  the  frequency  band,  with 
SNR  Gg/Gw  =  100  (20  dB).  The  signals  and  noise  are  band-pass  Gaussian 
random  processes.  The  sampling  rate  in  the  simulation  is  fg  =  10^  sam¬ 
ples/  s,  with  T  =  3  X  10^  samples,  so  the  time  interval  length  is  T  =  3  s. 

Figure  9a  displays  the  simulated  RMS  error  on  TDE  for  0.2  <  y*  <  1.0, 
along  with  the  corresponding  CRB  from  (25).  The  simulated  RMS  error  is 
based  on  100  runs,  and  the  TDE  is  estimated  from  the  location  of  the  maxi¬ 
mum  of  the  cross  correlation  of  the  sensor  signals.  The  threshold  coherence 
for  this  case  is  0.41,  from  (30)  and  (29).  Note  in  Figure  9a  that  the  simu¬ 
lated  RMS  error  on  TDE  diverges  sharply  from  the  CRB  very  near  to  the 
threshold  coherence  value  of  0.41,  illustrating  the  accuracy  of  the  analytical 
threshold  coherence  in  (30). 

Next  we  consider  TDE  with  three  different  signals: 

•  A  narrowband  signal  with  A/  =  2  Hz  centered  at  /o  =  40  Hz.  We 
refer  to  this  as  "1  harmonic"  (it  is  the  fundamental  frequency  of  the 
signals  defined  next). 
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Figure  9.  Comparison  of 
simulated  RMS  error  for  io 
TDE  wifh  CRBs  and 
threshold  coherence 
value,  (a)  Wideband 
signal  with  A/  =  30  Hz  g 
centered  at  /o  =  100  Hz.  io 
(b-d)  Narrowband 
signal  with  A/  =  2  Hz, 
fundamental  frequency 
/o  =  40  Hz  and  1,  2,  and 
5  harmonic 
componenfs, 
respectively. 


RMS  TDE  (Wideband):  =  100  Hz,  Af  =  30  Hz 


(b) 


RMS  TDE:  1  Harmonic,  L  =  40  Hz,  A  f  =  2  Hz 
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Coherence 


•  "2  harmonics"  at  40  and  80  Hz,  with  bandwidth  A/  =  2  Hz  at  each 
harmonic. 

•  "5  harmonics"  at  40,  80, 120, 160,  200  Hz,  with  bandwidth  A/  =  2  Hz 
at  each  harmonic. 

The  signal,  noise,  and  coherence  are  flat  over  each  frequency  band,  with 
SNR  Ga/Gw  =  100  (20  dB),  and  the  signals  and  noise  are  band-pass  Gaus¬ 
sian  random  processes.  The  sampling  rate  in  the  simulation  is  fa  =  10^ 
samples/ s,  with  T  =  2  x  10^  samples,  so  the  time  interval  length  is  T  =  2  s. 

Figures  9b-d  display  the  simulated  RMS  error  on  TDE  (based  on  1000  runs) 
for  coherence  values  0.7  <  7^  <  1.0.  As  in  the  previous  example,  the  TDE  is 
obtained  by  cross  correlation.  The  threshold  coherence  is  defined  only  for 
the  "1  harmonic"  signal,  and  the  threshold  coherence  value  is  «  1.  Figure  9b 
illustrates  the  divergence  of  the  simulated  RMS  error  from  the  CRB,  except 
at  7s  =  1. 

Figures  9c  and  d  display  the  results  for  the  signals  with  "2  harmonics"  and 
"5  harmonics."  The  additional  harmonics  enable  accurate  TDE  for  lower 
coherence  values,  but  we  cannot  use  (30)  to  compute  the  analytical  thresh¬ 
old  coherence  for  the  harmonic  signals.  The  "approximate  threshold  coher¬ 
ence"  values  indicated  in  Figures  9c  and  d  are  computed  as  follows.  For  the 
K  harmonics,  suppose  the  total  bandwidth  of  all  harmonics.  A/  =  K  -  2  Hz, 
is  centered  at  the  fundamental  frequency  /o  =  40  Hz.  The  approximate 
threshold  coherence  values,  0.95  in  Figure  9c  and  0.50  in  Figure  9d,  are 
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considerably  lower  than  the  actual  points  of  divergence  from  the  CRB.  Not 
surprisingly,  the  total  bandwidth  that  is  "spread"  across  the  harmonics  is 
less  useful  for  overcoming  signal  coherence  loss  than  an  equivalent  band¬ 
width  concentrated  at  the  fundamental  /o  =  40  Hz.  Narrowband  and  har¬ 
monic  signals  are  generally  difficult  for  TDE  due  to  ambiguous  peaks  in 
the  cross-correlation  function. 


3.3  TDE  With  H  >  2  Sensors 


We  can  extend  the  analysis  of  the  H  =  2  sensor  case  to  TDE  with  H  >  2 
sensors  following  the  approach  of  Weinstein  [18],  leading  to  the  conclusion 
that  pairwise  TDE  is  essentially  optimum  for  cases  of  interest  with  reason¬ 
able  signal  coherence  between  sensors.  By  pairwise  TDE  we  mean  that  one 
sensor,  say  H,  is  identified  as  the  reference,  and  only  the  H  —  1  time  differ¬ 
ences  Dih,  D2h,  ■  ■  ■ ,  Dh-i,h  are  estimated.  Under  the  conditions  described 
next,  these  H  —  1  estimates  are  nearly  as  accurate  for  source  localization  as 
forming  all  pairs  of  TDEs  Dg^  for  all  g  <  h.  Weinstein's  analysis  [18]  is  valid 
for  moving  as  well  as  nonmoving  sources. 

Extending  (27)  and  (28)  to  H  >  2  sensors,  let  us  assume  equal  Gg^hh/Gw  at 
all  sensors  h  =  1,. . .  ,H  and  equal  coherence  7s  between  all  sensor  pairs, 
so  that  the  SNR(|7s|)  in  (28)  is  equal  for  all  sensor  pairs.  Then  as  long  as 
H  ■  SNR(|7s|)  »  1,  improves  the  source  localization  variance  relative  to 
pairwise  processing  by  the  factor 


g(i  +  2-i7:) 


(31) 


Clearly  V  ^  1  as  7s  ^  1,  and  V  <  (3iE)/[2(l  -|-  H)]  <  1.5  for  jg  >  0.5. 
Therefore  the  potential  accuracy  gain  from  processing  all  sensor  pairs  is 
negligible  when  the  coherence  exceeds  the  threshold  values  that  are  typi¬ 
cally  required  for  TDE. 

This  result  suggests  strategies  with  moderate  communication  bandwidth 
that  potentially  achieve  nearly  optimum  localization  performance.  The  ref¬ 
erence  sensor,  iE,  sends  its  raw  data  to  all  other  sensors.  Those  sensors 
h  =  1, . . . ,  Tf  —  1,  locally  estimate  the  time  differences  Di^h,  . . . ,  Dh-i,h, 
and  these  estimates  are  passed  to  the  fusion  center  for  localization  pro¬ 
cessing  with  the  bearing  estimates  ^i, . . . ,  (j)H-  A  modified  scheme  with 
more  communication  bandwidth  but  more  centralized  processing  is  for  all 
H  sensors  to  communicate  their  data  to  the  fusion  center,  with  TDE  per¬ 
formed  at  the  fusion  center. 


3.4  TDE  With  Measured  Data 
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Eirst,  we  present  an  illustration  based  on  processing  the  measured  data  for 
the  source  in  Eigure  4a  that  was  discussed  in  section  2.2.2.  Eigure  10  shows 
results  of  cross-correlation  processing  of  the  data  for  a  2-s  segment  at  time 


Cross-Correlation:  Array  1,  Sensors  1,5 


Figure  10.  (a)  Cross 
correlation  of  signals  at 
arrays  1  and  3  in  Figure 
4a  for  source  at  time 
342  s.  (b)  Cross 
correlation  of  signals  at 
two  sensors  on  array  1, 
separated  by  8  ft. 


Lag  (s) 


Lag  (s) 


342  s.  Figure  10a  is  obtained  by  cross-correlating  the  signals  received  at  ar¬ 
rays  1  and  3,  for  which  the  coherence  is  appreciable  only  over  a  narrow 
band  near  39  Hz  (see  Figure  4c).  A  peak  in  the  cross  correlation  is  not  evi¬ 
dent,  which  is  expected  based  on  the  preceding  analysis,  since  nearly  perfect 
coherence  is  needed  for  narrowband  time  delay  estimation  in  this  scenario. 
Figure  10b  is  obtained  by  cross  correlating  the  signals  received  at  two  sen¬ 
sors  on  array  1,  where  the  coherence  is  large  over  a  wide  bandwidth  (see 
Figure  4d).  The  peak  is  clearly  evident  in  the  cross  correlation  in  Figure  10b. 

Next,  we  present  a  TDE  example  based  on  measured  data  from  Sanders 
Corporation  [19]  using  a  synthetically  generated,  nonmovmg,  wideband 
acoustic  source.  The  PSD  of  the  source  is  shown  in  Figure  11,  which  m- 


Figure  11.  (a)  Frequency 
spectrum  of  wideband 
signal,  (b)  Location  of 
nodes  where 
transmitters  and 
receivers  are  placed. 

(c)  PSD  and  coherence 
for  S5mthetically 
generated  wideband 
source  located  at  node  2 
and  measured  at  nodes 
0  and  1 .  (d)  Generalized 
cross  correlation  of  1-s 
segments  from  node  0 
with  a  10-s  segment 
from  node  1 .  The  data  is 
time-aligned  so  the  true 
peak  location  is  zero  lag. 


Frequency  (Hz) 
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dicates  that  the  source  bandwidth  is  about  50  Hz  with  center  frequency 
100  Hz.  With  reference  to  the  sensor  locations  in  Figure  lib,  the  source  is  at 
node  2  and  the  two  receiving  sensors  are  at  nodes  0  and  1.  The  source  and 
sensors  form  a  triangle,  with  dimensions  as  follows:  the  distance  from  the 
source  (node  2)  to  sensors  0  and  1  is  233  ft  and  329  ft,  respectively,  and  the 
distance  between  sensors  0  and  1  is  233  ft.  The  PSD  and  coherence  magni¬ 
tude  estimated  from  1-s  segments  of  data  measured  at  sensors  0  and  1  is 
shown  in  Figure  11.  Note  that  the  PSDs  of  the  sensor  signals  do  not  have 
their  maxima  at  100  Hz  due  to  the  acoustic  propagation.  However,  the  co¬ 
herence  magnitude  is  roughly  0.8  over  a  50-Hz  band  centered  at  100  Hz. 

Figure  12  shows  the  threshold  coherence  computed  with  (29)  and  (30)  for 
the  signal  in  Figure  11a  that  is  centered  at  a;o/27r  =  100  Hz  and  T  =  1-s 
observation  time.  For  bandwidth  At(;/27r  =  50  Hz,  the  threshold  coherence 
in  Figure  12  is  approximately  0.5.  The  actual  coherence  of  0.8  in  Figure  11c 
exceeds  the  threshold  value,  so  TDE  between  sensors  0  and  1  should  be 
feasible.  Figure  lid  shows  that  the  generalized  cross  correlation  has  its  peak 
at  zero  lag,  which  is  the  correct  location  because  the  sensor  data  is  time- 
aligned  before  processing.  This  example  shows  the  feasibility  of  TDE  with 
acoustic  signals  measured  at  widely  separated  sensors,  provided  that  the 
SNR,  fractional  bandwidth,  time-bandwidth  product,  and  coherence  meet 
the  required  thresholds. 


Figure  12.  Threshold 
coherence  for  S5mthetic 
wideband  source  with 
PSD  in  Figure  11a, 
which  has  center 
frequency  ujq  =  27rl00 
rad/s  and  observation 
time  is  T  =  1  s. 


(Oq  =  2  n  1 00  rad/sec,  T  =  1  s 
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Figure  13  contains  another  example  from  the  Sanders  data  using  the  sen¬ 
sor  locations  in  Figure  11b  and  the  wideband  source  with  spectrum  in  Fig¬ 
ure  11a.  In  this  example,  the  source  is  at  node  0,  and  the  receiving  sensors 
are  at  nodes  1  and  3.  Note  the  difference  in  the  PSD  shapes  in  Figure  13a, 
which  is  similar  to  our  observation  about  the  PSDs  in  Figure  11c.  The  signal 
coherence  between  nodes  1  and  3  is  shown  in  Figure  13b,  indicating  high 
coherence  over  an  appreciable  bandwidth.  The  cross  correlation  is  shown 
in  Figure  13,  and  the  peak  is  clearly  evident  at  the  correct  location. 

Figure  14  contains  a  final  example  using  the  Sanders  data  in  which  the 
source  is  at  node  0  and  measurements  are  recorded  at  nodes  1,  2,  and  3 
(see  Figure  11b  for  the  node  locations).  Differential  TDs  are  estimated  using 
the  signals  measured  at  nodes  1,  2,  and  3,  and  the  TDs  are  hyperbolically 
triangulated  to  estimate  the  location  of  the  source.  Figure  14a  shows  the  hy¬ 
perbolas  obtained  from  the  three  differential  TD  estimates,  and  Figure  14b 
shows  an  expanded  view  near  the  intersection  point.  The  triangulated  lo¬ 
cation  is  within  1  ft  of  the  true  source  location,  which  is  at  (—3.0)  ft. 

We  conclude  this  section  with  an  example  based  on  measured  data  in  an 
open  field.  The  source  is  an  Ml  tank  that  is  moving  at  a  range  of  approx¬ 
imately  140  m  from  a  set  of  H  =  3  sensor  arrays.  The  sensor  arrays  are 
labeled  A,  B,  and  C,  and  they  are  located  along  a  line  with  the  following 
separations: 


B 


15  m 


8  m 


C. 


Figure  13.  (a)  PSDs  at 
nodes  1  and  3  when 
transmitter  is  at  node  0. 

(b)  Coherence  between 
nodes  1  and  3. 

(c)  Generalized  cross 
correlation  between 
nodes  1  and  3. 
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Figure  14. 

(a)  Intersection  of 
hyperbolas  obtained 
from  differenfial  TDs 
esfimafed  af  nodes  1, 2, 
and  3  in  Figure  lib. 

(b)  Expanded  view  of 
part  (a)  near  the  point  of 
infersecfion. 


(b)  Hyperbolic  Triangulation  of  Time-Delays  from  Processed  Data 


The  source  is  moving  parallel  to  the  line  connecting  the  three  arrays.  Fig¬ 
ure  15  shows  the  PSDs  at  each  array  in  (a),  the  signal  coherence  between 
pairs  of  arrays  in  (b),  and  the  cross  correlations  in  (c).  The  PSDs  in  Fig¬ 
ure  15a  exhibit  strong  harmonics,  perhaps  due  to  tread  slap.  The  coherence 
in  Figure  15b  is  high  over  a  rather  large  bandwidth,  so  the  cross-correlation 
functions  in  Figure  15c  have  a  clear  peak  at  the  correct  location. 

The  examples  presented  in  this  section  have  demonstrated  success  and  fail¬ 
ure  of  time  delay  estimation,  depending  on  the  coherence  and  bandwidth 
of  the  measured  signals. 


Figure  15.  Results  for 
Ml  fank  moving  af 
range  140  m  from  fhree 
arrays  A,  B,  and  C 
separated  by  8, 15,  and 
23  m.  (a)  Power  spectral 
densities,  (b)  Signal 
coherence  between 
pairs  of  arrays,  (c)  Cross 
correlafions  between 
pairs  of  arrays. 
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4.  Subspace  Processing 


In  this  section,  we  begin  with  an  eigenanalysis  of  the  cross-spectral  density 
matrix  (14)  for  the  case  oi  H  =  2  arrays  containing  Ni  =  N2  =  N  sen¬ 
sors  operating  in  a  narrow  band  of  frequencies  centered  at  tco-  This  analysis 
leads  to  a  MUSIC -like  subspace  algorithm  for  source  localization  with  dis¬ 
tributed  arrays  and  partial  signal  coherence.  We  illustrate  with  computer 
simulations  and  measured  aeroacoustic  data  that  the  new  subspace  algo¬ 
rithm  is  limited  in  performance  due  to  source  location  ambiguities  that 
arise  from  the  large  separation  between  arrays  and  the  narrowband  signals. 
These  ambiguities  are  identical  to  those  that  we  studied  with  the  Ziv-Zakai 
bounds  for  TDE  in  section  3.1. 

4.1  Eigenanalysis  and  Algorithm 

To  simplify  notation,  we  will  let  a/,  represent  the  array  manifold  Uhiu^o) , 
represent  the  average  signal  power  at  array  h  in  the  frequency  band  of  in¬ 
terest,  is  the  average  noise  power,  and  7  is  the  coherence  7s,i2(<^o)-  We 
allow  7  to  be  complex-valued  with  I7I  <  1.  Then  (14)  can  be  expressed  as 
the  correlation  matrix 

t>(12)  /  0-f  \  aiaf  ai  a|^  exp(-j/UoT>i2)  Oj 

VL7Vi(T2  cr|  J  /  La2af  exp(-hjUoZli2)  a2a^ 

(32) 

where  ®  denotes  Kronecker  product,  o  denotes  elementwise  product,  and 
IvfxAf  is  an  A  X  A  matrix  of  I's.  The  superscript  (12)  in  (32)  indicates  that 
the  quantity  is  the  combined  correlation  matrix  of  all  sensors  in  arrays  1 
and  2.  The  correlation  matrices  for  the  individual  arrays  1  and  2  are 

=(T2a;,af  +  h  =  l,2-  (33) 

It  is  well  known  that  the  individual  array  correlation  matrices  RT)  and 
R(2)  in  (33)  have  a  rank-1  signal  subspace,  and  that  the  signal  subspace 
eigenvalue  is 

\f^=Nal  +  al,  h  =  l,2,  (34) 

with  corresponding  eigenvector 

h  =  l,2.  (35) 

The  remaining  eigenvalues  n  =  2, . . . ,  N  are  equal  to 

Now,  we  present  the  eigenanalysis  of  the  combined  correlation  matrix  RT^) 
in  (32).  The  derivation  of  these  results  is  similar  to  the  reasoning  in  the  ap¬ 
pendix  of  [20].  We  assume  that  the  arrays  are  labeled  so  that  af  >  cj|,  and 
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that  the  sensors  are  omnidirectional  so  that  the  array  manifold  elements  are 
of  the  form  [aLh\n  =  exp(— jcaor/i„).  The  rank  of  the  signal  subspace  compo¬ 
nent  of  is  at  most  two,  and  the  two  largest  eigenvalues  are 


=  y  (‘^i  +  •’I) 


1  +  ,  1-(1- 


\ 


Af  =  y  W  +  <^i) 


1  - 


1-(1- 


(T1  I  (T2 

(72  '  (^1  , 


£1  _|_  £2 
(72  '  (71  , 


The  corresponding  eigenvectors  are 


.(12) 


e  = 


1  I  £1  _  ^12 

2  \  0-2  CTl 


1  +  t  /  1  + 


2|7l 
£1_  £2 
0-2  CTl 


ai 


.(12) 


0O  — 


7*  exp(jwoT>i2)  a2 
-7  ai 


1  (  £1  _  £2 
2  V  CT2  cri 


1  +  t  /  1  + 


2|7l 
^1  ^2 
<72  cri 


+  cjI  (36) 


+  <7^  (37) 


exp(jwo-Di2)  a2 


provided  ui  /  <72-  Some  observations  are  as  follows. 

•  In  the  case  that  ui  =  <72  =  as,  (36)-(39)  reduce  to 


(38) 


,(39) 


(1  +  ItI  )  + 


= 


A 


^^2)  =  (1  -  IaI  )  + 


ai 

e(i2)  _ 

-(7/I7I)  ai 

.  (7VI7I)  exp(jwoT’i2)  a2  _ 

5  ^2 

_  exp^juioDu)  a2  _ 

(40) 

(41) 


Note  that  (41)  is  valid  as  long  as  7  /  0,  and  in  this  case  the  eigen¬ 
vectors  are  independent  of  the  coherence  magnitude.  In  the  case  of 
high  SNR  so  that  <C  1V<7^(1  —  I7I),  the  ratio  of  signal  eigenvalues 
is  approximately 


A 


(12) 


A 


(12) 


1  +  2 


I7I 


1  -  I7I 


(42) 


From  (42),  in  order  for  A^^  ’  >  IOA2  so  that  A^^^^  is  the  dominant 
eigenvalue,  the  coherence  must  satisfy  I7I  >  0.8181. 

In  the  case  of  perfect  coherence  I7I  =  1,  has  a  rank-1  signal  sub¬ 
space,  with  signal  eigenvalue  and  eigenvector 


+  12, CO J?)  Ar<  2  ,  2n  ,  2  (12,COH) 

A)  ’  '  =  +fj2)  +<7^,  e)  ’ 


<7i  ai 

[  <727*  exp(jwo-Di2)  a2 


(43) 
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The  expressions  in  (43)  are  the  natural  extension  of  (34)  and  (35)  ap¬ 
plied  to  the  "superarray"  of  all  sensors,  allowing  for  different  signal 
power  levels  at  the  two  arrays. 

•  In  the  case  of  partial  coherence  jyl  <  1,  has  a  rank-2  signal  sub- 

(12)  (12) 

space.  The  signal  eigenvalues  are  ordered  with  equality 

if  and  only  if  7  =  0  and  ai  =  a2-  In  the  incoherent  case  7  =  0,  the 
signal  eigenvectors  are 


(12, INC) 


ai 

(12, INC)  _ 

.  tdo  — 

0 

0 

*  Z 

.  ^2  _ 

(44) 


Thus,  in  the  incoherent  case,  the  signal  eigenvectors  (44)  of  are 
completely  determined  by  the  signal  eigenvectors  (35)  of  individual 
array  correlation  matrices  R*^^^  and  R^^). 


•  Note  that  the  span  of  the  signal  subspace  eigenvectors  in  (38) 
and  (39)  is  equal  to  the  span  of  the  incoherent  eigenvectors  in  (44). 
Therefore,  for  jyj  <  1,  the  signal  subspace  of  is  spanned  by  the 
signal  subspace  eigenvectors  of  the  individual  array  correlation  matri¬ 
ces  R(^)  and  R(^).  It  follows  that  R^^^)  and  the  block  diagonal  ma¬ 


trix 


have  identical  signal  and  noise  subspaces.  Thus  a 


R(i)  0 

0  R(2) 

subspace  algorithm  that  exploits  joint  processing  of  the  two  arrays 
must  exploit  the  structure  of  the  dominant  signal  subspace  eigenvec¬ 
tor  in  (38). 


The  following  subspace  procedure  may  be  used  to  estimate  the  source  lo¬ 
cation  {xs,  Us)  using  the  correlation  matrices  R^^^,  R^^^,  and 

1.  Form  the  correlation  matrices  R^^^  and  R^^^  using  the  observations 
from  the  individual  arrays.  Estimate  the  noise  power  from  the 
noise  eigenvalues  of  R^^^  and  R^^^.  Then  estimate  the  signal  power 
at  each  array  using  the  largest  eigenvalue  of  R^^^  and  R^^^  with  (34): 

^  (Af ^  >  h  =  l,2.  (45) 


2.  Form  the  joint  correlation  matrix  and  estimate  the  coherence 

magnitude  I7I  using  the  two  largest  eigenvalues  A^^^\  A^^^^  of  R(^^) 
with  (36)  and  (37): 


ItI 


1  - 


1  - 


.(12)  ,(12) 
N{af+a'^) 


1/2 


(46) 


as  long  as  ui  /  cr2.  Note  that  if  ai  =  a2,  then  it  is  not  necessary  to  Es¬ 
timate  the  coherence  magnitude  (see  (41)  and  step  3).  We  will  assume 
that  7  is  real-valued  and  positive,  which  is  motivated  by  physical  con¬ 
siderations  as  discussed  earlier  and  in  Wilson  [7, 8]. 
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3.  Use  the  estimates  for  ai,o'2,  and  I7I  in  the  signal  eigenvector  expres- 

(12) 

sion  for  ^  in  (38)  to  create  a  "search  vector"  g: 

/5i  I  (47^ 


g{xs,ys)  = 

1  (  £1  _  £2  ) 

2  \<J2  cri  y 

ai 

* 

L  7 

exp(yu;oUi2)  a2 

where  ai,  a2,  and  D12  are  functions  of  the  source  location  {xs,  Vs)-  (If 
cji  =  CJ2,  then  the  search  vector  is  based  on  e)  ^  in  (41)  instead.) 

4.  Let  en‘^\  n  =  1, . . .  ,2N  denote  the  eigenvectors  of  ordered 

(12)  (12) 

such  that  the  corresponding  eigenvalues  follow  A)  ^  >  •••  >  ^2N  ■ 
Then  form  the  matrix 


V  = 


e(i2) 

^2N 


(48) 


A  MUSIC-type  spectrum  may  then  be  defined  as 

Pc{Xs,ys)  =  g/ryv^g’ 

and  the  value  of  {xs,ys)  that  maximizes  (49)  is  an  estimate  of  the 
source  location. 

5.  Wax  and  Kailath  [21]  proposed  an  "incoherent"  MUSIC-type  source 
location  estimator  that  ignores  coherence  between  the  arrays.  The 
form  of  the  Wax  and  Kailath  estimator  can  be  expressed  as  follows: 

Pi(xs,ys)  = - u - u — 1  (50) 

af  V(i)v(i)^ai  -h  af  V(2)v(2)^a2 

where  . . .  ,e^^]  contains  the  noise  eigenvectors  of 

This  analysis  extends  to  the  case  of  H  arrays  with  a  single  source. 
If  the  signals  are  partially  coherent  at  all  arrays,  then  the  dimension 
of  the  signal  subspace  of  the  combined  correlation  matrix  is  H.  If  K 
signals  are  present,  then  the  combined  correlation  matrix  has  signal 
subspace  with  rank  KH. 


Examples  of  Subspace  Processing 

We  present  examples  of  the  "partially  coherent"  MUSIC  algorithm  (49)  and 
the  "incoherent"  MUSIC  algorithm  (50)  with  computer  simulations  and 
measured  data.  First,  we  simulate  a  scenario  with  H  =  2  arrays.  The  in¬ 
dividual  arrays  are  identical  and  contain  Ni  =  N2  =  N  =  7  sensors.  Each 
array  is  circular  and  has  4-ft  radius,  with  six  sensors  equally  spaced  around 
the  perimeter  and  one  sensor  in  the  center.  Narrowband  processing  cen¬ 
tered  at  40  Hz  is  assumed,  with  an  SNR  of  40  dB /sensor  at  array  1  and 
34  dB  at  array  2,  i.e.,  =  10^,  Gs,22{^)/Gw{u})  =  2500.  The 

coherence  is  set  to  two  values,  7^72  =  0.9  and  0.999.  The  arrays  are  located 


at  coordinates  (xi,  yi)  =  (0, 150),  (x2, 2/2)  =  (0, 0),  and  one  source  is  located 
at  {xs,ys)  =  (200, 300),  where  the  units  are  meters.  Source  location  estima¬ 
tion  is  performed  using  M  =  30  snapshots,  where  each  snapshots  contains 
an  estimate  of  the  measured  complex  amplitude  at  40  Hz  at  each  sensor. 

Results  are  presented  in  Figure  16,  where  part  a  shows  an  approximate 
CRB  as  a  function  of  signal  coherence.*  Increasing  coherence  leads  to  mod¬ 
est  improvement  in  the  localization  potential.  Figures  16b  and  16c  contain 
representative  spectra  of  incoherent  and  partially  coherent  MUSIC,  respec¬ 
tively,  for  coherence  7,;, 12  =  0.9.  Note  that  the  incoherent  MUSIC  spectrum 
has  a  rather  broad  but  unique  maximum.  The  partially  coherent  MUSIC 
spectrum  has  a  sharp  ridge  in  one  direction,  with  multiple  peaks  along  the 
ridge.  The  multiple  peaks  arise  due  to  the  ambiguities  in  source  location 
caused  by  the  narrowband  signals  and  large  separation  between  arrays. 


Figure  16.  (a)  CRBs  on 
localization  accuracy, 
(b)  Incoherent  MUSIC 
and  (c)  partially 
coherent  MUSIC 
spectra. 


Incoherent  MUSIC  (WAX/KAILATH)  Partially  Coherent  Music 
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Table  1  contains  results  of  Monte  Carlo  simulations  with  100  runs  for  both 
coherence  values,  7^72  =  0.9  and  0.999.  The  incoherent  MUSIC  performs 
close  to  the  approximate  CRB  in  Figure  16a  corresponding  to  7^72  =  0.  Par¬ 
tially  coherent  MUSIC  performs  worse  than  incoherent  MUSIC  because  of 
the  ambiguities:  partially  coherent  MUSIC  sometimes  chooses  the  incorrect 
peak  along  the  ridge  in  the  spectrum  shown  in  Figure  16c.  Note  that  par- 

*The  CRB  is  approximate  because  the  mapping  of  the  scaling  factor  in  (16)  to  the 
number  of  snapshofs  M  is  approximate. 
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Table  1.  RMS  error  simulation  results  for  incoherent  and  partially  coherent  MUSIC. 


Coherence 

Approx 

CRB 

Incoherent  MUSIC 

Partially  Coherent  MUSIC 

Xs 

Vs 

Xs 

Vs 

Xs 

Vs 

0.9 

1.17 

1.10 

1.36 

1.10 

2.14 

1.71 

0.999 

0.81 

0.87 

1.41 

1.16 

1.96 

1.61 

0.0 

1.50 

1.22 

tially  coherent  MUSIC  is  more  accurate  for  js  =  0.999  than  for  js  =  0.9, 
suggesting  that  the  ambiguities  are  less  severe  when  the  signal  coherence 
is  larger. 

We  have  applied  incoherent  MUSIC  and  partially  coherent  MUSIC  to  the 
measured  data  scenario  discussed  in  section  2.2.2  and  in  Figures  4  and  10. 
The  data  is  processed  in  a  narrow  band  around  39  Hz  for  the  2-s  interval 
between  345  and  347  s.  Both  algorithms  produce  identical  estimates  for  the 
source  location  in  this  case,  but  the  spectra  for  partially  coherent  MUSIC  are 
shown  in  Figure  17.  The  top  panel  clearly  shows  the  ambiguities  that  limit 
the  performance  with  narrowband  signals.  The  contour  plot  in  the  bottom 
panel  also  shows  the  periodic  structure  of  the  ambiguities. 


Figure  17.  Partially 
coherent  MUSIC 
spectrum  for  measured 
dafa. 


(Og  =  2  71  50  rad/s,  T  =  2  s 
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5.  Moving  Sources 


One  of  our  objectives  in  this  report  is  to  quantify  scenarios  in  which  jointly 
processing  data  from  widely  spaced  sensors  has  the  potential  for  improved 
source  localization  accuracy,  compared  with  incoherent  triangulation/ 
tracking  of  bearing  estimates.  We  established  in  section  2  that  the  poten¬ 
tial  for  improved  accuracy  depends  directly  on  TDE  between  the  sensors. 
Then  we  showed  in  section  3  that  TDE  with  partially  coherent  signals  is  fea¬ 
sible  only  with  an  increased  time-bandwidth  product  of  the  sensor  signals. 
This  leads  to  a  constraint  on  the  minimum  observation  time,  T,  in  passive 
applications  where  the  signal  bandwidth  is  fixed.  If  the  source  is  moving, 
then  approximating  it  as  nonmoving  becomes  poorer  as  T  increases,  so 
modelling  the  source  motion  becomes  more  important. 

Approximate  bounds  are  known  [22,  23]  that  specify  conditions  of  validity 
for  nonmoving  and  moving  source  models.  Eet  us  consider  H  =  2  sensors 
with  Doppler  values  02  >  oi  (see  (64)  for  the  definitions  of  oi,  02)-  If  /max 
(Hz)  is  the  maximum  signal  frequency  that  is  processed,  then  TDE  estima¬ 
tion  accuracy  is  not  seriously  affected  by  ignoring  source  motion,  as  long 
as  the  time  interval  T  satisfies 


1 


fn 


(51) 


Taking  typical  parameters  for  ground  vehicles  in  aeroacoustics,  let  us  con¬ 
sider  a  vehicle  moving  at  5%  the  speed  of  sound  (15  m/ s),  with  radial  mo¬ 
tion  that  is  in  opposite  directions  at  the  two  sensors.  Then  02/01  —  1  ~  0.1 
and  (51)  becomes  T  <C  10//max-  For  /max  =  100  Hz,  the  requirement  is 
T  <C  0.1  s,  which,  according  to  the  analysis  in  section  2,  yields  insufficient 
time-bandwidth  product  for  partially  coherent  signals  that  are  typically  en¬ 
countered.  Thus  motion  modeling  and  Doppler  compensation  are  critical, 
even  for  aeroacoustic  sources  that  move  more  slowly  than  in  this  example. 


5.1  Data  Model 


In  this  section,  we  extend  the  nonmoving  source  model  from  section  2  us¬ 
ing  first-order  motion  models  (see  (52),  (53),  (67)).  The  first-order  motion 
models  are  simple  and  accurate  over  larger  time  intervals  T  compared 
with  the  nonmoving  source  model.  However,  accurate  modeling  of  more 
complex  trajectories  over  longer  time  intervals  requires  higher-order  poly¬ 
nomial  models,  with  added  complexity.  The  source  position  trajectory  is 
modeled  as  a  straight  line  with  constant  velocity  over  an  interval  of  length 
T, 

Xs{t)  =  Xsfl  +  Xs-  {t-to),  to<t<to  +  T  (52) 

ys{t)  =  ysfl  +  Vs-  {t-to),  (53) 
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so  Xs,  its  are  the  velocity  components.  The  source  trajectory  parameter  vec¬ 
tor  is 

&  =  [xsfl,Xs,ysfl,ysY' ,  (54) 

and  the  (time-varying)  propagation  time  from  the  source  to  the  sensors  on 
array  h  follows  from  (1)  and  (2): 


Thit) 

Thn{t) 


=  -  [iXs{t)  -  Xhf  +  {ys{t)  -  yhf] 
cc 


Xs{t)  -Xh  ys{t)  -yh  . 

^Xhn  H - T-7T^ - ^Vhr, 


dh{t) 


dh{t) 


(55) 


=  [{cos  (j)h{t))  Ax hn  +  {sin  (j)h{t)) Ay hn]  ■  (56) 

The  bearing  and  bearing  rate  are  related  to  the  source  motion  parameters 
0  as 


4’h{t) 

^h{t) 


tan 


-1 


ys{t)  -  Vh 

Xs{t)  -  Xh 


=  tan 


-1 


ijs  cos  (t)h{t)  -  Xs  sin  (t)h{t) 


dh{t) 


ys,o  +  ys-  {t-  to)  -  yh 
Xsfi  +  Xs-  {t-  to)  -  Xh 


(57) 

(58) 


The  radial  velocity  of  the  source  with  respect  to  array  h  is 


Vr,h{t)  =  XsCOS(j)h{t)  +  ys^dn(t>h{t)- 


(59) 


We  can  insert  (52)  and  (53)  into  (55)  to  obtain  the  following  approximation 
for  the  propagation  time  to  array  h: 


=  Th{to) 


^  2  •  cos  (/>fe(to)  -  Xs  -  {t-  to)  2  •  sin  <j)h{to)  ■  ys  ■  {t  -  to)  ] 


dh{to) 


dh{to) 


Th{to)  +  •  {t  -  to), 


(60) 

(61) 


where  dh{to)  and  Vr,h{to)  are  the  source  distance  and  radial  velocity  at  the 
start  of  time  interval  t  =  to-  The  approximation  (61)  is  valid  as  long  as  the 
total  motion  during  the  time  interval  T  is  much  less  than  the  range,  i.e., 

\2xsT\  <C  dh{to)  and  |2y^r|  <C  dh{to). 

Next,  we  use  the  approximation  (61)  and  model  the  received  signal  at  the 
reference  sensor  on  array  h  as 


Sh{t-Th{t))  =  Sh 


1  -  V  -  r,(fo)  + 


Sh{aht-  Th{to)  -6  )  ,  to<t<to  +  T, 


(62) 

(63) 


where 


^r,h{io) 


=  1 - [xs  COS  (ph{to)  -6  ijs  sin  (ph{to)] 


ah  =  l 


c 


c 


(64) 


is  the  Doppler  compression  and 

Tfc(to)  =  =  -  [{xs,o  -  Xhf  +  iVsfi  -  Vhf]  (65) 

c  c 

is  the  propagation  delay  at  the  initial  time  t  =  to.  Without  loss  of  generality, 
we  set  to  =  0,  so  the  received  signal  at  sensor  n  on  array  h  is 

Sh  (ah  t  -  r/i(0)  -  Thn{t))  ,  (66) 

which  is  the  extension  of  the  signal  component  of  (8)  to  the  moving  source 
case.  Note  from  (2)  that  Thn{t)  depends  on  the  source  location  only  through 
the  time-varying  bearing  (j)h{t),  which  we  approximate  with  a  first-order 
model 

(ph{t)  ^  (ph{to)  +  ^{to)  -{t-to),  to  <t  <to  +  T.  (67) 

For  a  single  array  h,  the  Doppler  compression  ah  and  time  delay  Th{to) 
have  negligible  effect  on  estimation  of  the  intra-array  delays  Thn{t),  since 
ah  and  Th{to)  are  identical  for  each  n  =  1, . . . ,  Nh-  Thus,  each  array  can  be 
processed  separately  to  estimate  the  bearings  (t>i{to), . .  • ,  (t>H{to)  and  bear¬ 
ing  rates  0i(fo);  •  •  •  >  ^H{to),  and  these  can  be  "triangulated"  via  (57)  and 
(58)  to  estimate  the  source  motion  parameters  0  in  (54).  An  algorithm  [24] 
for  estimating  0/i(io)  and  ^h{to)  is  described  in  section  5.2. 

Let  us  consider  the  signals  received  at  the  reference  sensors  at  each  array, 
so  Thn{t)  =  0  in  (66): 

t  -  Ti{to)  ^1  -  ^  ^  _  Tnito)  ■  (68) 

Our  modeling  assumptions  imply  that  each  signal 
Sh  (^1  —  f  _  Th{to)  is  a  wide-sense  stationary  Gaussian  ran¬ 

dom  process.  However,  for  two  arrays  g,  h  with  unequal  Doppler 
Xr,g{to)  /  Xr^h{to),  the  signals  at  arrays  g,  h  are  not  jointly  wide-sense 
stationary  [22,  25],  complicating  the  analytical  description  and  the  CRB 
performance  analysis.  The  jointly  nonstationary  sensor  signals  generally 
are  not  characterized  by  a  cross-spectral  density  matrix,  TDE  with  jointly 
nonstationary  signals  as  in  (68)  was  derived  by  Knapp  and  Carter  [22].  The 
CRB  analysis  was  rigorously  justified  by  Schultheiss  and  Weinstein  [26], 
and  they  extended  the  results  to  CRBs  on  differential  Doppler.  A  clever 
transformation  is  used  so  that  the  jointly  nonstationary  signals  in  (68)  are 
locally  modeled  by  a  CSD  of  the  form  (14),  and  Schultheiss  and  Weinstein 
show  [25]  that  the  representation  is  accurate  for  CRB  analysis. 

Formulating  the  results  of  Schultheiss  and  Weinstein  [25]  for  the  case  of 
partially  coherent  signals,*  we  make  the  following  observations  for  H  =  2 
arrays.  The  results  are  valid  for  large  observation  time  (T  much  larger  than 
the  coherence  time  of  the  signals  and  noise).  The  TDE  is  D12  =  n  (fo)  — ''‘2(^0) 
and  the  differential  Doppler  is  Avi2  =  Vr,i(to)  —  Vr,2{to)- 

*The  signal  coherence  between  the  signals  at  arrays  g  and  h  in  (68)  is  defined  assuming 
perfect  compensation  of  the  Doppler  compression  ag ,  an,  thus  yielding  the  definition  in  (7). 
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•  Estimation  of  TDE  and  differential  Doppler  are  decoupled,  so  the 
CRB  on  Di2  is  given  by  (24),  which  is  identical  to  the  nonmoving 
source  case. 

•  The  threshold  coherence  analysis  for  TDE  in  (30)  and  Eigure  8  extends 
to  the  moving  source  case.  In  the  best  case  that  Doppler  effects  are 
perfectly  estimated  and  compensated,  the  TDE  problem  that  remains 
is  identical  to  the  nonmoving  source  case.  Doppler  estimation  is  less 
demanding  in  terms  of  time-bandwidth  product  compared  with  TDE. 
Indeed,  Doppler  estimation  is  possible  with  sinusoidal  signals  that 
have  negligible  bandwidth  [25]. 

The  CRB  on  differential  Doppler  [25],  modified  for  partially  coherent 
signals,  is 


CRB(Az;i2) 


247r  /  c  \  2 

~T~  \t) 


-1 


SNR  (I75, 12(0;)!)  d(jj 


(69) 


Note  that  (69)  is  a  scalar  multiple  of  the  CRB  on  TDE  in  (24).  The 
CRB  on  differential  Doppler  may  be  achievable  in  scenarios  where 
the  time-bandwidth  product  is  insufficient  for  TDE. 

Interestingly,  differential  Doppler  provides  sufficient  information  for 
source  localization,  even  without  TDE,  as  long  as  five  or  more  sen¬ 
sors  are  available  [26].  Thus  the  source  motion  may  be  exploited  in 
scenarios  where  TDE  is  not  feasible,  such  as  narrowband  signals. 

•  We  discussed  TDE  with  H  >  2  sensors  in  section  2,  concluding  that 
pairwise  processing  of  TDEs  Dih,  . . . ,  Dh-i,h  with  a  reference  sen¬ 
sor  H  is  nearly  optimum  for  scenarios  of  interest  (see  (31)).  A  similar 
result  holds  for  differential  Doppler  estimation  [26],  where  pairwise 
estimation  of  Avih,  ■  ■  ■ ,  Avh-i,h  is  nearly  as  accurate  as  estimation 
of  all  pairs  Avgh,  as  long  as  H  •  SNRdy^l)  »  1. 


5.2  An  Algorithm 

The  parameters  that  can  be  directly  estimated  from  the  sensor  data  are  the 
bearings  ipiito),  ■  ■  ■,  bearing  rates  0i(fo);  •  •  •  >  pairwise  time 

differences  Dm  =  ri(fo)  -  Tuito),  ■  ■  ■,  Dh-i,h  =  TH-i{to)  -  rnito),  and 
differential  Doppler  Avm  =  Vr,i{to) -Vr,H{to),  ■■■,  Avh-i,h  =  Vr,H-i{to)  - 
Vr,H{to)-  Equations  (57),  (58),  (64),  and  (65)  define  the  nonlinear  relations 
that  "triangulate"  these  parameters  and  relate  them  to  the  source  motion 
parameters  0  =  [xsfl,Xs,ys,o,ys]'^- 

A  distributed  processing  algorithm  is  outlined  next,  and  parts  of  the  algo¬ 
rithm  are  illustrated  with  measured  aeroacoustic  data  in  section  6. 

1.  Use  the  local  polynomial  approximation  (EPA)  beamformer  [24]  at 
each  array  to  estimate  the  bearings  and  bearing  rates.  The  EPA  beam- 
former  in  Katkovnik  and  Gershman  [24]  is  formulated  for  narrow- 
band  processing,  and  it  is  a  generalization  of  the  classical  beamformer 
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to  moving  sources.  We  extend  it  in  a  straightforward  way  to  wide¬ 
band  signals  by  incoherently  averaging  the  LPA  beampatterns  at  dif¬ 
ferent  frequencies. 

2.  Solve  (57)  and  (58)  to  obtain  initial  estimates  of  the  source  motion  pa¬ 
rameters  0.  These  estimates  correspond  to  incoherent  triangulation 
of  the  bearings  and  bearing  rates  from  individual  arrays. 

3.  Estimate  the  Doppler  compression  factors  ui , . . . ,  an,  compensate  for 
Doppler,  and  test  whether  the  signals  at  distinct  arrays  have  sufficient 
coherence,  fractional  bandwidth,  and  time-bandwidth  product  to  en¬ 
able  TDE  between  arrays  (see  section  2  for  the  conditions). 

4.  If  the  conditions  are  not  met,  then  incoherent  triangulation  of  the  bear¬ 
ings  and  bearing  rates  is  nearly  optimum,  and  further  joint  processing 
is  not  informative. 

5.  If  the  conditions  are  met,  then  identify  a  reference  array  H  (the  ar¬ 
ray  with  maximum  SNR)  and  estimate  the  time  differences  Dih,  . . ., 
Dh-i,h  arid  differential  Dopplers  Avih,  . . .,  Avh-i,h-  The  maximum 
likelihood  solution  involves  wideband  ambiguity  function  search  over 
Doppler  and  TDE  [22],  but  computationally  simpler  alternatives  have 
been  investigated  [26]. 

6.  A  suboptimum  procedure  is  to  avoid  the  joint  Doppler  and  TDE  es¬ 
timation  in  the  preceding  step,  and  instead  use  the  initial  Doppler 
estimates  from  steps  1  and  2  and  perform  TDE  after  approximate 
Doppler  compensation.  With  this  approach,  triangulation  of  the  TDEs 
via  (65)  will  improve  the  estimates  of  arid  ys,o  only  (and  not  the 
source  velocity  Xs,ys)- 

7.  If  multiple  sources  are  present,  then  the  EPA  beamformer  in  step  1 
may  be  used  to  separate  the  source  signals  at  each  array  prior  to 
Doppler /TDE  estimation. 

The  LPA  beamformer  in  steps  1  and  7  is  illustrated  in  section  6  for  a  two- 
source  scenario  based  on  measured  aeroacoustic  data.  Examples  of  TDE 
(step  5)  with  measured  data  were  presented  in  section  3.4. 
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6.  Example  with  Two  Moving  Sources 


We  present  an  application  of  the  LPA  beamformer  [24]  to  measured  aeroa- 
coustic  data  with  two  ground  vehicles:  an  Ml  tank  and  an  M3  tracked 
armored  vehicle.  The  vehicle  trajectories  over  a  10-s  segment  are  shown 
in  Figure  18a.  Two  arrays,  labeled  8B  and  8C,  are  separated  by  about  23 
m.  Each  array  is  circular  with  N  =  7  sensors,  4-ft  radius,  and  six  sen¬ 
sors  equally  spaced  around  the  perimeter  with  one  sensor  in  the  center. 
We  present  results  based  on  processing  the  data  at  array  8B,  for  which  the 
bearing  and  bearing  rate  of  the  sources  are  shown  in  Figure  18b.  The  bear¬ 
ing  of  the  M3  varies  by  more  than  30°  over  the  10-s  time  interval.  The  range 
of  the  M3  is  closer  at  about  100  m,  with  the  Ml  range  at  approximately  200 
m. 

The  data  from  array  8B  is  processed  over  a  wide  bandwidth  from  30  to 
150  Hz.  The  beampattern  for  a  classical  beamformer  based  on  a  nonmoving 
source  model  is  shown  in  Figure  18c.  The  beampattern  is  the  incoherent 
sum  of  narrowband  beampattern  over  the  30-  to  150-Hz  frequency  band. 
The  peak  of  the  beampattern  is  located  at  approximately  the  mean  bearing 
of  the  stronger  source  (M3)  over  the  10  s  interval. 

The  beampattern  of  the  FPA  beamformer  is  shown  in  Figure  19a.  The  FPA 
beampattern  exploits  a  first-order  model  for  time-varying  bearing,  as  in 
(67).  The  FPA  beampattern  is  two-dimensional,  with  axes  of  initial  bearing 
and  bearing  rate  <^(to)-  The  FPA  beamformer  in  Katkovnik  and  Ger- 
shman  [24]  is  formulated  for  narrowband  processing,  so  we  extend  to  the 
wideband  case  by  incoherently  adding  the  narrowband  FPA  beampatterns 
over  the  frequency  range.  The  peak  of  the  beampattern  in  Figure  19a  is  close 
to  the  true  values  of  (pito),  p{to)  for  the  M3  that  are  shown  in  Figure  18b. 
The  location  of  the  weaker  Ml  source  is  not  evident  in  the  FPA  beampattern 
in  Figure  19a,  so  we  subtract  an  estimate  of  the  stronger  source  from  the 
data.  The  subtraction  is  performed  based  on  the  bearing  and  bearing  rate 
estimates  from  Figure  19a.  The  subtraction  is  coherent  over  the  process¬ 
ing  bandwidth  and  includes  the  time-varying  bearing.  The  FPA  beampat¬ 
tern  after  subtraction  is  shown  in  Figure  19b,  which  indicates  the  bearing 
and  bearing  rate  of  the  weaker  Ml  source.  This  example  illustrates  the  gain 
in  resolution  that  is  achieved  by  exploiting  the  source  motion  in  a  beam- 
former. 
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Bearing  B 


Figure  18.  (a)  Array 
locations  and  trajectory 
of  Ml  and  M3  vehicles 
over  a  10-s  time 
interval,  (b)  Bearing  and 
bearing  rate  of  sources 
with  respect  to  array  8B. 
(c)  Beampattern  of 
classical  beamformer 
based  on  nonmoving 
source  model. 
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Figure  19. 

(a)  Beampattern  of  LPA 
beamformer,  showing 
the  location  of  the 
stronger  source  (M3). 

(b)  Beampattern  of  LPA 
beamformer  after 
subtracting  an  estimate 
of  the  strong  source 
signal  from  the  data, 
showing  the  location  of 
the  weaker  source  (Ml). 
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7.  Summary  and  Concluding  Remarks 


We  have  presented  in  this  report  an  analysis  of  source  localization  with  sen¬ 
sors  arranged  in  an  "array  of  arrays"  configuration.  We  have  paid  partic¬ 
ular  attention  to  aeroacoustical  localization  of  ground  vehicles,  where  the 
signals  measured  at  widely  separated  sensors  are  not  perfectly  coherent 
due  to  random  propagation  effects.  We  analyzed  an  algorithm  that  com¬ 
bines  with  pairwise  TDE  between  separate  arrays.  We  showed  that  this 
scheme  incorporates  distributed  processing  and  data  compression  so  that 
the  communication  bandwidth  with  a  fusion  center  is  reduced,  with  little 
loss  in  localization  accuracy.  The  algorithm  was  presented  first  for  a  single, 
nonmoving  source,  and  then  we  extended  to  the  case  of  moving  sources. 
The  algorithm  for  moving  sources  involves  estimation  of  bearing-rate  and 
differential  Doppler  between  sensors.  We  provided  an  analysis  based  on 
Ziv-Zakai  bounds  that  quantifies  the  requirements  on  SNR,  signal  band¬ 
width,  signal  coherence,  and  observation  time  so  that  joint  (coherent)  pro¬ 
cessing  of  widely  spaced  sensor  data  provides  improved  localization  ac¬ 
curacy.  We  presented  computer  simulations  and  results  from  processing 
measured  data  to  illustrate  and  support  the  theoretical  developments  in 
the  report. 

We  are  continuing  our  research  on  aeroacoustic  source  localization  and 
tracking  in  several  ways.  In  particular,  we  are  investigating  algorithms  that 
exploit  the  "sum  of  harmonics"  structure  that  characterizes  the  acoustical 
signature  of  many  ground  vehicles.  We  are  also  investigating  joint  algo¬ 
rithms  for  localization,  tracking,  and  classification  that  share  information 
to  improve  performance. 
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